suppressPackageStartupMessages({
library(tidyverse)
library(pheatmap)
library(here)
library(ggiraph)
library(ggrepel)
})
Differential Expression (DTE & DGE) Heatmaps
top_tx_df <- tx_df %>%
mutate(is_sig = FDR_DEG < 0.05, abs_logFC = abs(logFC_DEG)) %>%
arrange(desc(is_sig), desc(abs_logFC))
top100_tx <- top_tx_df %>% slice_head(n = 100) %>% pull(row_id)
top50_tx <- head(top100_tx, 50)
top_gene_df <- gene_df %>%
mutate(is_sig = FDR_DEG < 0.05, abs_logFC = abs(logFC)) %>%
arrange(desc(is_sig), desc(abs_logFC))
top100_genes <- top_gene_df %>% slice_head(n = 100) %>% pull(row_id)
top50_genes <- head(top100_genes, 50)
plot_tx_heatmap <- function(tx_ids, title) {
mat <- cpm_matrix[rownames(cpm_matrix) %in% tx_ids, , drop = FALSE]
mat <- log2(mat + 1)
mat <- mat[rowSums(mat) > 0, , drop = FALSE]
if (nrow(mat) == 0) return(NULL)
# Use the modified label map with PB accessions and asterisks
rownames(mat) <- tx_label_map[rownames(mat)]
annotation_df <- data.frame(
sample_name = colnames(mat),
group = ifelse(grepl("WT", colnames(mat)), "WT", "Q157R")
) %>% column_to_rownames("sample_name")
annotation_colors <- list(group = c(Q157R = "#336B87", WT = "#E99787"))
col_order <- c(grep("Q157R", colnames(mat), value = TRUE),
grep("WT", colnames(mat), value = TRUE))
mat <- mat[, col_order]
annotation_df <- annotation_df[col_order, , drop = FALSE]
n_rows <- nrow(mat)
plot_height <- min(10, 5 + n_rows * 0.15)
fontsize_row <- if (n_rows <= 30) 10 else if (n_rows <= 60) 8 else 6
pheatmap(mat, scale = "row", cluster_rows = TRUE, cluster_cols = FALSE,
annotation_col = annotation_df, annotation_colors = annotation_colors,
show_colnames = TRUE, show_rownames = TRUE, main = title,
fontsize_row = fontsize_row, fontsize_col = 10, height = plot_height)
}
plot_gene_heatmap <- function(gene_ids, title) {
mat <- gene_cpm[rownames(gene_cpm) %in% gene_ids, , drop = FALSE]
mat <- log2(mat + 1)
mat <- mat[rowSums(mat) > 0, , drop = FALSE]
if (nrow(mat) == 0) return(NULL)
rownames(mat) <- gene_label_map[rownames(mat)]
annotation_df <- data.frame(
sample_name = colnames(mat),
group = ifelse(grepl("WT", colnames(mat)), "WT", "Q157R")
) %>% column_to_rownames("sample_name")
annotation_colors <- list(group = c(Q157R = "#336B87", WT = "#E99787"))
col_order <- c(grep("Q157R", colnames(mat), value = TRUE),
grep("WT", colnames(mat), value = TRUE))
mat <- mat[, col_order]
annotation_df <- annotation_df[col_order, , drop = FALSE]
n_rows <- nrow(mat)
plot_height <- min(10, 5 + n_rows * 0.15)
fontsize_row <- if (n_rows <= 30) 10 else if (n_rows <= 60) 8 else 6
pheatmap(mat, scale = "row", cluster_rows = TRUE, cluster_cols = FALSE,
annotation_col = annotation_df, annotation_colors = annotation_colors,
show_colnames = TRUE, show_rownames = TRUE, main = title,
fontsize_row = fontsize_row, fontsize_col = 10, height = plot_height)
}
plot_tx_heatmap(top50_tx, "Top 50 Differentially Expressed Transcripts")
plot_tx_heatmap(top100_tx, "Top 100 Differentially Expressed Transcripts")
plot_gene_heatmap(top50_genes, "Top 50 Differentially Expressed Genes")
plot_gene_heatmap(top100_genes, "Top 100 Differentially Expressed Genes")
Differential Transcript Usage (DTU) Heatmaps
top_tx_dtu_df <- tx_dtu_df %>%
arrange(adj.p.value_DTU, desc(abs(lf_DTU)))
top100_tx_dtu <- top_tx_dtu_df %>% slice_head(n = 100) %>% pull(row_id)
top50_tx_dtu <- head(top100_tx_dtu, 50)
top_gene_dtu_df <- gene_dtu_df %>%
arrange(adj_p.value_DTU, desc(abs(lf_DTU)))
top100_gene_dtu <- top_gene_dtu_df %>% slice_head(n = 100) %>% pull(row_id)
top50_gene_dtu <- head(top100_gene_dtu, 50)
plot_dtu_heatmap <- function(mat, row_label_map, row_ids, title) {
mat <- mat[rownames(mat) %in% row_ids, , drop = FALSE]
mat <- mat[rowSums(mat) > 0, , drop = FALSE]
mat <- mat[apply(mat, 1, function(x) sd(x, na.rm = TRUE) > 0), , drop = FALSE]
if (nrow(mat) == 0) return(NULL)
rownames(mat) <- row_label_map[rownames(mat)]
annotation_df <- data.frame(
sample_name = colnames(mat),
group = ifelse(grepl("WT", colnames(mat)), "WT", "Q157R")
) %>% column_to_rownames("sample_name")
annotation_colors <- list(group = c(Q157R = "#336B87", WT = "#E99787"))
col_order <- c(grep("Q157R", colnames(mat), value = TRUE),
grep("WT", colnames(mat), value = TRUE))
mat <- mat[, col_order]
annotation_df <- annotation_df[col_order, , drop = FALSE]
n_rows <- nrow(mat)
plot_height <- min(10, 5 + n_rows * 0.15)
fontsize_row <- if (n_rows <= 30) 10 else if (n_rows <= 60) 8 else 6
pheatmap(mat, scale = "row", cluster_rows = TRUE, cluster_cols = FALSE,
annotation_col = annotation_df, annotation_colors = annotation_colors,
show_colnames = TRUE, show_rownames = TRUE, main = title,
fontsize_row = fontsize_row, fontsize_col = 10, height = plot_height)
}
plot_dtu_heatmap(frac_matrix, frac_label_map, top100_tx_dtu, "Top 100 DTU Transcripts")
plot_dtu_heatmap(frac_matrix, frac_label_map, top50_tx_dtu, "Top 50 DTU Transcripts")
plot_dtu_heatmap(gene_frac_matrix, gene_frac_label_map, top100_gene_dtu, "Top 100 DTU Genes")
plot_dtu_heatmap(gene_frac_matrix, gene_frac_label_map, top50_gene_dtu, "Top 50 DTU Genes")